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1  Introduction 


In  ([8])  and  ([30]),  we  introduce  a  new  approach  called  the  Method  of  Auxiliary  Map¬ 
ping,  to  deal  with  domain  singularity  and  interface  singularity,  which  arise  in  such 
elliptic  boundary  value  problems  as  steady  state  heat  transfer.  In  this  paper,  this 
approach  is  extended  for  plane  elasticity  problems  containing  singularities. 

There  are  three  versions  of  the  finite  element  method:  the  h-version,  the  p-version, 
and  the  h-p  version.  The  h- version  ([13], [35])  is  the  standard  one,  where  the  degree  p 
of  the  elements  is  fixed,  usually  at  a  low  level,  typically  with  p  =  1,2,  or  3  and  the 
accuracy  is  achieved  by  properly  refining  the  mesh.  The  p-version  ([10],  [12],  [36]), 
in  contrast,  fixes  the  mesh  and  achieves  better  accuracy  by  increasing  the  degree  p 
of  the  elements  uniformly  or  selectively.  The  k  .'ron  ([4]-[5],  [11],  [16]-[21])  is  a 
combination  of  both.  In  this  paper,  we  are  con»>  vith  the  p-version  of  the  finite 
element  method. 


In  the  theory  and  practice  of  the  finite  element  method,  much  work  has  been  done 
to  design  special  approaches  to  deal  with  elasticity  problems  containing  singularities 
([25], [32], [33], [38]).  Singularities  occur  when  the  solution  domain  has  cornerp  abrupt 
changes  in  boundary  data,  or  consists  of  two  or  more  materials.  These  siugularities 
are  called  a  comer  singularity  ([9], [15], [27], [34]),  a  boundary  data  singularity([29],[3‘7]), 
and  an  interface  singularity([23],[26],[29],[31]),  respectively. 

In  an  effort  to  provide  accurate  and  economical  solutions,  many  different  approaches 
to  deal  with  singularity  in  elasticity  problems  have  been  attempted  over  the  years  Basi¬ 
cally  there  are  three  ways  the  problem  is  approached:  mesh  refinement([7],[ll],[16],(17j, 
[18],  [19],  [20], [35]);  use  of  special  singular  elements([l],[2],[22]);  and  use  of  (nonlocal) 
special  singular  functions([24]).  Expanding  the  trial  space  by  adding  special  singular 
(local  or  global)  functions  which  mimic  the  singularities  can  lead  to  a  more  accurate  so¬ 
lution,  but  more  problems  will  be  generated,  especially  in  computer  coding.  Moreover, 
one  must  know  the  structure  of  the  eigenvalues  corresponding  to  the  singular  points 
in  order  to  choose  proper  singular  functions.  The  most  popular  approach  is  the  mesh 
refinement,  but  its  success  depends  on  a  proper  choice  of  mesh  and  it  also  requires 
longer  computing  time.  Moreover,  when  the  singularity  is  very  strong,  as  in  Example 
5.II,  this  approach  cannot  give  any  acceptable  results. 

In  this  paper,  the  Method  of  Auxiliary  Mapping  introduced  in  ([8])  and  ([30])  will  be 


modified  so  that  it  can  efficiently  handle  the  singularities  which  arise  in  plane  elasticity 
problems.  It  will  be  shown  that  this  new  approach  yields  far  better  results  for  elasticity 
problems  containing  singularities  than  do  conventional  approaches  at  virtually  no  extra 
cost.  Moreover,  this  method  gives  a  reasonable  solution  for  those  elasticity  problems 
which  even  can  not  be  solved  by  using  the  h-p  version  of  the  Finite  Element  Method. 

This  paper  is  organized  as  follows:  The  notations  and  the  model  problems  to  work 
with  are  described  in  §2.  In  §3,  the  structure  of  the  comer  singularity  and  basic  lemmas 
are  introduced.  In  §4,  the  Method  of  Auxiliary  Mapping  is  explained  in  the  context  of 
plane  elasticity,  and  the  improvement  of  error  bounds  by  Method  of  Auxiliary  Mapping 


□n 


is  presented.  Various  numerical  results  to  demonstrate  effectiveness  of  our  method  are 
given  in  §5.  These  include  an  especially  remarkable  success  in  Examples  5.1  and  5.II 
and  MAM’s  handling  of  interface  singularities  caused  by  an  abrupt  change  in  material 
properties.  Finally,  the  concluding  remarks  are  given  in  §6. 


2  Preliminary 

2.1  The  Notation 

For  £71^  a,  polygonal  domain  with  boundary  dCl,  we  let  L^iVl)  = 
k  >  0  integer,  denote  the  usual  Sobolev  spaces.  For  u  €  we  denote  by  ||u|ljfe,n 

2tnd  |u|ik,n,  the  usual  norm  and  seminorm,  respectively. 

In  elasticity,  the  state  variable  is  the  displacement  vector  denoted  by 
{u}  =  {ux(a!,y),Uy(x,y)}^  and  the  flux  is  the  stress  tensor  denoted  by  {<7^“!}  = 

Let  strciin  tensor.  Then  the  strain- 

displacement  and  the  stress-strain  relations  axe  given  by 

{£<“)}  =  [Cl{u)  and  {<r(">}  =  |f;j{e<“))  (1) 


respectively,  where  [jD]  is  the  differential  operator  matrix: 


i  dy  dx  \ 


and  [E\  =  [J5,j],l  <  t,j  <  3,  is  the  symmetric  positive  definite  matrix  of  material 
constants.  For  an  isotropic  elastic  body,  in  the  case  of  plane  stress. 


[E] 


E 

1-1/2 


1 

u  1 
0  0 


0 

0 

1-1/ 

2 


where  E  is  the  module  of  elasticity  and  i/(0  <v<  1/2)  is  Poisson’s  ratio. 
The  equilibrium  equations  of  elasticity  are 

y)  -I-  {/}(i, y)  =  0,  (x,  y)  €  n. 


(2) 


where  {/}  =  {/*(x,  y),  fy{x,  y)}^  is  the  vector  of  internal  sources  representing  the  booy 
force  per  unit  area. 
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2.2  The  Model  Problem 


Introducing  the  relations  (1)  into  (2),  the  equilibrium  equations  can  be  expressed  in 
terms  of  the  displacement  vector  {u}.  Consider  the  following  system  of  the  partial 
differential  equations  in  terms  of  the  displacement  vector: 

(Dn£][D]{u}(i,y)+ {/}(x,y)  =  0,  {x,y)en  (3) 

subject  to  the  following  boundary  conditions: 

[;v]{ffW}(^)  =  =  >er',  (4) 

{“}(“)  =  {“}(«)=  s€r“,  (5) 

where  U  =  dfl,  {n*,  riy}^  is  a  unit  vector  normal  to  the  boundary  dil  of  the 
domain  fl  and 

1^1  =  [  0  n. 

Let  =  {{to}  =  {tOx,tOy}  €  :  {to}  =  0  on  F^}.  Then,  as  usual, 

the  variational  form  of  (3)-(5)  is  as  follows:  find  the  vector  {«}  such  that  Ux,Uy  € 
{u}  =  {u}  on  F^  and 

Biiuh  {o})  =  for  all  {v}  €  Hh{n),  (6) 

where 

8(W.{''})  =  [m{v)f[E]{\D]{u})dxdy,  (7) 

Ja 

^({u})  =  J^{vy{f}dxdy  +  j)^{vy{f}d3.  (8) 

By  the  strain  energy  of  the  displacement  vector  {u}  we  mean  W({u})  =  (l/2)B({u},  {u}). 

The  finite  element  approximation  of  the  solution  of  (6)  is  to  construct  approxima¬ 
tions  of  each  component  of  the  vector  {u}.  We  denote  the  basis  functions  defined  on  ft 
by  '9i{x,y),i  =  1,2,  The  components  of  the  displacement  vector  in  term  of  basis 
functions  are  of  the  form: 

n 

i=l 

n 

»=i 

where  a,(t  =  1,2,  ...,2n)  are  called  the  amplitudes  of  the  basis  functions  Let 

{*.}  =  I  =  (9) 
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Then  {«}  caji  be  written  as: 


{u}  = 

«=i 

Moreover  we  have  the  following: 

Lemma  2.1  The  bilinear  form  B({u},{v})  on  an  element  e  becomes 

^  (V*j)  ./{o}  =  {«.,0)’'(.nrf{tx)  =  {»„0)’', 

jf(V»,f  ^  (V«i)  i/W  =  {0,*.f  aniiW  =  {«j,0)’'. 


3  The  Corner  Singularities 

The  accuracy  of  the  finite  element  approximation  depends  on  the  regularity  of  the  true 
solution  ([6]).  In  the  presence  of  singularity  the  solution  of  (3)-(5)  has  a  low  regularity. 
In  this  section  the  structure  of  the  singularity  due  to  nonsmoothness  of  the  domain 
will  be  investigated.  For  this  purpose  the  equations  of  elasticity  (3)  will  be  localised 
by  restricting  (3)  on  a  neighborhood  of  a  singulaxity  of  the  domain  fi. 


3.1  The  behavior  of  a  solution  in  the  vicinity  of  nonsmooth 
boundary 

Let  us  consider  the  equations  of  elasticity  (3)  in  the  vicinity  of  a  corner  shown  in 
Fig.  3.1.  When  the  body  force  is  neglected,  the  equations  of  elasticity  in  the  polar 
coordinates  system  can  be  written  as 

/>  «  ^  ^  ldu0  Ur,  Id  (due  IdUr  , 

« + +  ;  W  +  7>  -  ;  W  +  7>  =  “  <“> 

« + 2");  ar + ;  as- + 7  J + '‘a;<  ar  - ;  W + 7> = « 


where 


E  ,  _  vE 
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which  characterize  materials.  Here  E  is  the  module  of  elasticity  amd  v  is  Poisson’s 
ratio. 

The  radial  and  tangential  stresses  on  the  wedge  surfaces  are 


dur  \  dug  ,  .1  due  Ur . 

.due  ,  Idur  ue. 

;  IF '  7>- 


(14) 

(15) 


Let  us  consider  the  solutions  of  (11)  and  (12)  in  the  following  form: 

X,  =  7/(«), 


(16) 

(17) 


Substitution  of  these  forms  into  (11)  and  (12)  gives  a  system  of  ordinary  differential 
equations  in  /  and  g.  One  can  see  that  the  solution  of  this  system  has  the  form 

f  =  Ai  cos  {(1  +  A)5]  +  A2  sin  {(1  +  A)5] 

+i43Cos((l  -  A)5]  +  A4sin[(l  -  A)d] ,  (18) 

g  —  A2  cos  ((1  +  A)0]  —  Ai  sin  [(1  +  A)^] 

+»//l4  cos  [(1  —  A)0]  —  fiAz  sin  [(1  -  A)0] ,  (19) 

where  =  (3  +  A  —  4i/)/(3  —  A  —  4i/).  Then  displacements  and  stresses  in  the  vicinity 
of  the  comer  will  be  expressed  as 


=  j4i  cos  [(Id-  A)5]  -1-  i42  sin 

[(1  -h  X)9] 

-|-A3COs[(1  — 

A)0]  4-  A4  sin  [(1  —  A)0] , 

(20) 

r~^ue 

=  A2  cos  [(1  +  A)0]  —  Ai  sin 

1(1 + m 

d-»;444Cos[(l  • 

-  A)0]  —  i/Aasin  [(1  -  A)0] , 

(21) 

fi  r  (Tee 

=  — 2Ai4icos[(l 

+  A)^]  —  2Ay42  sin  [(1  -t-  A)^] 

-(1  +  A)(1- 

»/)y43COs[(l 

-m 

-(1  -h  A)(l  - 

i7)A|sin  [(1 

-  m , 

(22) 

H  r  CTrB 

=  — 2A44i  sin  [(1 

d"  A)0]  d"  2AA2  cos  [(Id"  A)0] 

-(1  -  A)(l  - 

ij)A3sin[(l 

-m 

-(1-A)(1- 

17)A4COs[(1 

-  A)^] . 

(23) 

From  now  on,  ’’the  wedge  angle  a”  stands  for  "the  wedge  angle  20”.  Let  us  label 
the  boundary  conditions  along  the  boundaries  0  =  a  and  9  =  —a  of  a  wedge-shaped 
region  bounded  by  radii  6  =  ±a,  as  shown  in  Fig.  3.1,  as  follows: 

BCi  :  uj  =  u,  =  0, 

BC2  •  (^eS  —  ‘^rS  —  0, 

BC3  :  ue  =  Tro  =  0. 
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In  order  to  determine  the  eigenvalues  A  and  the  constants  Ai,A2,  A3,  A4,  for  various 
boundary  conditions  on  the  boundary  surfaces  6  —  ±q,  these  boundary  conditions  are 
applied  to  the  equations  (20)-(23).  Then  because  of  the  given  boundary  conditions,  we 
obtaun  the  following  trigonometric  equations  to  determine  the  eigenvalues  A  (see,  [25] 
and  [32]  for  details): 


sin(2AQ) 

sin(2Aa) 

sin(2Aa) 

sin^(2Aa) 

sin(4Aa) 

sin(4Ao) 


±Asin(2Qf)/C 

±Asin(2a) 

±  sin(2Q) 

(1  +  A)V4A- A2sin2(2a)/C 
— Asin(4Q) 

A  sin(4Q!)/C’ 


if  BC\-BC\  is  imposed, 
if  BC2-BC2  is  imposed, 
if  BC3-BC3  is  imposed, 
if  BC1-BC2  is  imposed, 
if  BCi-BC^  is  imposed, 
if  BC1-BC3  is  imposed. 


where  C  =  3  —  A/v  for  plane  strain  and  C  =  (3  —  +  t')  for  the  plane  stress  and 

BCi-BCj  means  that  BCi  is  on  one  side  and  BCj  is  on  the  other  side  of  the  boundary 
surfaces,  9  =  ±q. 

These  trigonometric  equations  could  have  complex  roots  as  well  as  real  roots.  For 
various  wedge  angles  a  and  for  various  boundary  conditions,  min{i?e(A)}  are  computed 
in  (§3.8  of  [32]).  In  particular,  if  the  condition  BC3  is  imposed  on  both  of  0  =  ±a, 
0  <  a  <  90®,  then  the  smallest  real  eigenvalue  A  (min{i?e(A)})  can  be  arbitrary  small 
as  a  goes  to  90®.  Fig.  3.2  shows  the  smallest  real  root  A  of  the  trigonometric  equation; 
sin2Aa  =  sin  2a,  for  75®  <  a  <  90®.  In  such  cases,  the  singularities  are  too  strong  to 
obtain  any  reasonable  approximations  by  standard  numerical  approaches.  These  cases 
will  be  elaborated  in  Example  5. II. 


Fig.  3.1.  The  vicinity  of  a  singularity  with  wedge  angle  a. 
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Fig.  3.2.  The  smallest  eigenvalue  versus  the  wedge  angle  a  when  the  normal 
displacement  and  tangential  traction  are  zero  along  the  boundary  surfaces,  9  =  ±a. 


3.2  The  Computation  of  the  Bilinear  form  on  the  Regions 
transformed  by  Auxiliary  Mappings 

Now  we  consider  an  auxiliary  mapping  which  will  be  used  to  describe  our  new  approach 
to  deal  with  elasticity  problems  containing  singularities.  Let  z  =  x-i-iy  and  C  =  ^  +  *»/• 
Let 

S  =  {(r,0);  0<r<R,a<9<b},  (24) 

5*  =  {(r*,r):  0<r^  <R^f^,al^<9*  <bll3}  (25) 

be  two  circular  sectors  in  the  z-plane  and  the  ^-plane  respectively.  Suppose  :  S*  —*  S 

is  the  conformal  mapping  defined  by 

Z  =  /(C)  =  C^  (26) 


and  let  0  be  the  inverse  function  of  /,  then  the  determinants  of  their  Jacobians  are 


|J(/)|  =  and  IJW)!  =  (27) 


8 


respectively.  We  denote  the  shifted  function  onto  5*  of  a  function  /  :  S  — »  R  by  the 
conformal  mapping  hy  f  =  /  o 

The  following  lemma  was  proved  in  [30]  and  it  will  play  a  key  role  in  the  Method 
of  Auxiliary  Mapping. 

Lemma  3.1.  For  u,  v  €  H^{S),  we  have 


On  ai2 
oji  aj2 


Vvdxdy=  [  (Vuf 
Js* 


9i2 
921  922 


Vvdidri 


(28) 


where 

'  <  =  (1  - 

9ii  =  ^11  i  +  022  sin*  t  —  (021  +  O12)  sin  t  cos  t 
'  9i2  =  (oii  —  022)  sin  t  cost  —  021  sin*  <  + 012  cos*  t 

921  =  (oii  —  022)  sin  t  cos  t  —  012  sin*  t  +  021  cos*  t 

922  =  Oil  sin*  t  +  022 cos*  t  +  (012  +  021)  sin t  cost, 

and  represents  the  polar  coordinates  of  points  in  5*.  For  v  6  R^iS)  and  f  € 

H^{S),  we  have 

ff{x,y)vix,y)dxdy  =  f  +  V^)^~'^f{(,r))H(,‘n)d^dr].  (29) 

Js  Js* 


Remark  3.2.  (1)  From  Lemma  2.1,  the  bilinear  form  5({^j},  {^j})  for  the  basis 
vector  functions  has  the  same  form  as  the  left  hand  side  of  the  integrals  in  Lemma  3.1. 
Thus,  the  local  stiffness  matrices  on  S  by  the  basis  vector  functions  (9)-(10)  are  the 
same  as  the  local  stiffness  matrices  on  5*  by  the  mapped  basis  vector,  {'^i}  = 

(2)  Constraining  the  stiffness  matrices  on  the  elements  in  S  due  to  the  non-homogeneous 
essential  boundary  conditions  {tXx»«y}^  along  parts  of  the  boundary  dS  of  S  has  the 
same  effect  as  constraining  the  stiffness  matrices  obtained  by  the  right  hamd  side  of 
Lemma  3.1  by  using  the  boundary  conditions  {uxoy?^,  UyO(^}^  along  the  corresponding 
boundaries  of  S*. 


4  The  Method  of  Auxiliary  Mapping:  A  New 
Approach  to  Deal  with  Singulsirities 

4.1  Description  and  Implementation  of  the  Method 

Overall  structure  of  our  method  is  as  follows.  From  the  arguments  given  in  section  2.1, 
if  no  body  forces  are  present,  then  in  a  neighborhood  Sp  of  a  singular  point  P,  each 
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component  of  the  displacement  vector  can  be  written  as 


Q(P) 

KjFj(r)Gjie)  +  WQiP^ir,  6).  (30) 

i=i 

Here  (r,  0)  are  the  polar  coordinates  with  respect  to  the  singular  point  P.  The  func¬ 
tion  Gj{6)  is  analytic  up  to  the  boundary  of  Sp  and  Fj{r)  =  Re{r^^3  log'’(r))  or 
log'’(r)),  where  p  =  0  except  for  some  special  angles.  The  eigenvalues  \pj 
are  in  general  complex  numbers  with  positive  real  parts  and  Re{Xpj)  <  Re{Xpj+i). 
Fj{r),Gj{0)  depend  on  the  interior  (wedge)  angle  a;  Kj  are  stress  intensity  factors. 
wq^p)  is  smoother  than  the  first  term  on  the  right  hand  side  of  (30).  Q{P)  is  a  posi¬ 
tive  integer  and  Re(XpQ^f.fj  <  1.  Let  Xpj  be  the  smallest  real  number  of  Re(Xpj)  and 
suppose  Ap]  <  1,  If  we  let  ^  =  1/Ap]  and  (r*,5*)  denotes  the  polar  coordinate  system 
on  Sp  centered  at  P*,  then 

Q(P) 

np{r\0-)  =  E  KiF^(ir-r)Gji^e*)  +  wcnp){ir^f,m  (31) 

i=i 

and  /3Apj  >  1  for  j  =  1,...,(3(P).  Therefore  «p(r*,0*)  =  (up  o  v?^)(r*,fl'')  is  in 
/f”‘(5>),m  >  2. 

Thus,  on  Sp,  the  standard  finite  element  method  could  yield  a  good  approximation 
of  Up  in  H^(Sp).  Since  our  auxiliary  mappings  are  conforming  and  the  mapping  sizes 
^  are  assumed  to  be  >  1,  the  P*-norm  is  preserved  under  the  transformation  by  the 
auxiliary  mapping.  Thus,  approximating  up  in  H^{Sp)  by  the  standard  finite  element 
basis  functions  defined  on  Sp  has  the  same  effect  as  approximating  up  on  5p  in  H^{Sp) 
by  using  singular  basis  functions  on  Sp  constructed  through  the  auxiliary  mapping 

However,  the  novelty  of  our  method  lies  in  never  constructing  such  singular  basis 
functions. 

We  now  describe  our  method,  the  Method  of  Auxiliary  Mapping  (MAM).  Suppose 
the  exact  solution  {u}  has  singularities  at  Pi,  Pj?  •••»  Pm-  In  what  follows,  u  denotes  the 
displacements  Ux  and  Uy.  In  this  case,  our  method  goes  as  follows: 

Step  1:  Determination  of  the  Singular  Regions 

At  each  singular  point  P,,  construct  a  neighborhood  of  the  singular  point  P,,  a 
sector  5,  centered  at  P,.  Namely, 


Sq  =  {(r,^)  :0<r<  P,}nn, 

where  (r,  0)  are  the  polar  coordinates  at  P,.  Our  method  is  not  sensitive  to  the  size 
of  the  radius  Rg,  provided  P,  is  chosen  small  enough  so  that  5,  is  a  circular  sector  in 
n  and  any  two  different  neighborhoods  of  two  singularities  are  disjoint.  Rg  is  usually 
selected  to  be  <  1. 

Step  2:  Selection  of  Auxiliary  Mappings 
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Suppose  is  min{i?e(A,j)},  where  A„  are  the  eigenvalues  of  the  singularity  at 
P,.  Then  the  mapping  sizes  of  auxiliary  mappings  are  selected  as  follows: 


(32) 


Now  the  auxiliary  mapping  :  S*  — ♦  5,  is  defined  by  z  =  (p^*(C)  =  ^  conformal 

mapping  from  the  C-pl^e  to  the  2-plane. 

Step  3:  Triangulation  of  f) 

For  each  5,  ,  generate  a  curvilinear  triangulation  7^  of  5,  as  shown  in  Fig.  4.1. 
Then  construct  a  triangulation  T  on  fl  such  that  T(s,  =  7^.  Let  T*  be  the  image  of 
Tg  imder  (see,  Fig.  4.1). 

For  e*  €  7^*  (  c  €  To  =  T\  U7^),  $e«  (  $e)  is  the  usual  elemental  mapping  from  the 
standard  element  fl,t  (which  is  either  the  reference  triangle  or  rectangle  depending  on 
whether  c*  (c)  is  a  triangular  or  a  rectangular  element)  onto  curvilinear  elements  e*(e) 
respectively.  Since  we  allow  circular  arcs  as  sides  of  elements,  the  elemental  mappings 
could  be  of  the  blending  type([14])  as  those  in  chapter  6  of  ([36])  and  satisfy  the  usual 
technical  conditions([4],[20])  that  lead  to  conforming  finite  elements. 

Step  4;  Computation  of  Stiffness  Matrix  and  Load  Vector. 

In  computing  local  stiffness  matrices  and  local  load  vectors, 


•  Use  the  standard  elemental  mapping  for  the  elements  e  in  the  nonsingular 

region  flo  =  5,. 

•  Use  the  standard  elemental  mappings  $e*  for  the  elements  e  in  the  singular 

regions  in  other  words,  local  stiffness  matrices  juid  load  vectors  on  the 

element  e  in  the  singular  region  are  replaced  by  those  computed  on  the  elements 
c,  =  by  using  the  right  hand  sides  of  equations  (28)  and  (29). 


Let  denote  the  special  elemental  mapping  from  Qst  onto  e  €  7^  defined  by 
o  $e. .  We  will  call  this  special  elemental  mapping  the  singular  elemental 

mapping. 


Remark  4.1.  (1)  If  is  used  as  the  element^  mapping  on  the  element  e  in 
a  singular  region  5,  then  the  basis  functions  constructed  through  will  mimic  the 
original  singularity  on  5,. 

(2)  Suppose  ei  €  To  ,  e2  €  U7^  and  7  =  ei  n  ej  =  {(ro,0)  :  a  <  9  <  b}.  Then  the 
conformal  mapping  (y?^*)"*  is  linear  on  the  clc«ed  interval  [a,  6);  therefore,  the  basis 
functions  constructed  by  using  the  usual  elemental  mapping  for  e  C  flo  ^d  the 
singular  elemental  mapping  for  e  C  US,  are  continuous  along  their  common  edges. 
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Let  be  the  standard  shape  functions  on  and  Mj  =  A^  o($f 

Then  =  Mj.  Hence,  from  Lemma  2.1  ,  we  have 

|j[  =  Jl^Mng>MVAr;fd(  (33) 

Jlrndz  =  (34) 

Instead  of  computing  the  left-hand  sides  of  (33)-(34)  involving  singular  shape  functions, 
we  compute  the  right-hand  sides  of  the  equations  for  the  local  stiffness  matrix  and  load 
vector  on  the  elements  c  in  the  singular  regions.  Let  us  note  that  from  Lemma  3.1 
the  coefficient  qij  are  not  singular.  Moreover,  If  is  chosen  to  be  an  integer  >  2,  the 
integrand  of  the  left  hand  side  of  (34)  is  as  smooth  as  /. 

Thus,  the  computer  implementation  of  our  method  is  quite  simple  since  any  existing 
finite  element  code  can  be  used  for  the  computation  of  the  right-hand  sides  without  any 
alterations.  Indeed,  in  MAM,  unlike  other  singular  function  approaches,  the  banded 
structure  of  the  resulting  stiffness  matrix  is  not  lost  and  no  severe  problem  with  ill 
conditioning  will  occur. 


Fig.  4.1.  A  singular  neighborhood  5,  of  a  singular  point  P,  and  its  mapped 
domain  S*  under  the  mapping  The  scheme  on  5*  and  the  corresponding  mesh 
Tq  on  Sq. 
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Let  Vp  =  {{uj}  =  {wi,W2)  €  :  to,|e  o  $e{w>i|e  o  ^f)  is  a  polynomial  of  degree 

p  on  Qrt  for  all  elements  c  in  flo  (c  in  U^i5,)},  where  !),<  is  the  standard  triangle, 
T,  or  the  standard  rectangle,  Q,  according  to  whether  e  is  a  triangular  element  or  a 
rectangvilar  element.  Then  the  p-version  of  the  Finite  Element  Method  in  our 
context  is  as  follows:  Find  an  element  {«p}  6  Vp  such  that 

=  :F({u}),  for  all  H  €  Kp.  (35) 

The  dimension  of  Vp  will  be  denoted  by  Np  and  will  be  called  the  degree  of  free- 
dom(DOF).  Let  us  note  that  in  the  p-version  of  the  finite  element  method  the  trian¬ 
gulation  of  n  is  fixed  and  only  the  degree  p  of  the  basis  polynomials  is  increased. 

If  {tXe*}  is  the  solution  of  (6)  then 

||{up}  -  {u„}l|£  =  rnin  -  {ue*}||£;,  (36) 

MeVp 

where  ||{to}|||;  =  ^B({iu},  {tu})  is  the  energy  norm. 

A 

4.2  The  Rate  of  Convergence. 


Let  flo  and  5,  be  the  same  as  those  in  the  preceding  argument.  Suppose  Ue*|no  € 
ff*^(f^o),  «e*|s,  €  If‘'^(Sg),  Vo  >2,  Vg<2,  I  <q<  M,  and  if  =  min{ile(A„)}  < 

1,  then  our  method  with  auxiliary  mapping  =  1/A^i^  will  greatly  reduce  the  in¬ 
tensity  of  the  singularity  at  P,.  Therefore,  Ue®  =  o  and  v*  >  2, 

which  is  larger  than  i/,. 

By  using  the  inequality  (see  [30]  for  detail), 

||u||i,s  <  ||w||i,s*, 


the  following  theorem  was  proved  in  ([30]). 

Theorem  3.1.  Suppose  is  the  finite  element  solution,  on  a  quasi  uni¬ 

form  mesh,  obtained  by  employing  the  method  of  auxiliary  mapping  with  the  auxiliary 
mapping  9^’  on  each  singular  region  Sq  in  the  framework  of  the  p-version  of  the  finite 
element  method.  Then  we  have 


^  ||«ei||wB,no  ,  ||«ex||i/;.5; 


(37) 


where  Np  is  the  degree  of  freedom  and,  for  each  q,  0  <  q  <  M,  Cq  is  independent 
ofNp. 
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4.3  Computation  of  Strains  and  Stress 

The  strains  are  computed  from  the  strain-displacement  relationships  and  the  stresses 
are  computed  directly  from  the  strain-stress  relationships.  However,  in  MAM,  those 
computations  on  the  singular  regions  are  different  from  the  standard  approach. 


(4.3A).  Suppose  ^  is  the  inverse  mapping  of  the  auxiliary  mapping  :  S*  —*  S, 
defined  by  (26),  and  J{<^)  =  d(xuX2)/d{^i,(2)  denotes  the  Jacobian  matrix  of  <^. 
Then  [./(^)  o  <^j  •  =  I  and  hence 

[J(Vp)r=J(V»)o/.  (38) 

Moreover, 


J(V>)  =  J{tp)o<fi^ 


cos(^-^)^ 

8in(^-^)^ 


—  sin( 
cos( 


1-/9 


/9 


)0 

)9 


-  «n(i  -  /9)r  ■ 

j3^  ’  sin(l  -  cos(l  -  ' 


(39) 


Let  us  recall  the  singular  elemental  mapping  for  an  element  e  is  defined  by 

=  y)  0  e, 

where  (^  :  e*  — »  e  is  the  auxiliary  mapping  defined  by  z  =  and  ^  :  E  -*  e*  is  the 
standard  elemental  mapping.  Now  we  have 

[j(i>^)Y'  =  (J(v =«)!-' 

=  J([y>o$]”')o((^o$) 

=  [j($~^  0(^"^)j  o((^o$) 

=  [(*^(^"*)0V’~*)*«/(v>"^)]  o(v’o^) 

=  {(([A^)ro^-*)ov’-*)- 

((•^(V’)]”^  o  V»~')}  o{(po^)  by  (38) 

and  hence 

=  (Wv)]"  o  if  ■  (imf)-'-  (40) 
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Suppose  P  —  (arijXj)  is  a  point  in  an  element  e  C  and  let  $  ;  >  e  be  the 

elemental  mapping.  Then  the  strains  at  P  are 

■  1  0 1  r  0  0  ■ 

=  0  0  V,Ui(P)+  0  1  V^U2{P) 

.  0  1  J  1 1  0 . 

■  1  0  1 

=  0  0 

0  1  J  /=! 

‘0  0]  „ 

+  0  1  [J($r]'  (P)Sa/+„V(J\r/(P), 

1  0  J  /=! 

wh«e  t(P)  = />.  V.  =  ( A,  l-f  and  V,  =  . 

Suppose  e  is  an  element  in  the  singular  region  on  which  MAM  is  applied,  then 
is  replaced  by  and  it  follows  from  (39)  and  (40)  that 


=  (ij(9)i-'o«)’'-([j(0)r')’' 

_  1/  .Xi-^  cos(l  -  /9)5*  -  8in(l  -  i9)0* 
^  sin(i  -  cos(i  -  p)e* 


,  (41) 


where  J?  is  the  (z,j)  component  of  ([J($)]^)~*. 

It  should  be  noted  that  in  (41),  the  exponent  (1  —  /?)  <  0  since  the  mapping  size  P 
is  larger  than  1.  Thus,  vinlike  the  standard  finite  element  solution,  from  (41),  we  can 
see  that  the  strains  at  the  singular  points,  calciilated  from  the  finite  element  solution 
obtained  by  using  MAM,  are  infinity.  This  is  because  our  solution  by  MAM  resembles 
the  exact  solution  near  the  singularities. 

(4.3B)  The  stresses  are  computed  from  the  stress-strjun  relation.  Usually  the 
principal  stress  and  the  equivalent  stress  are  of  engineering  interest.  The  eigenvalues 
of  the  following  matrix 

<r\\  <ri2 

<7ii  Oi2 

are  called  the  principal  stresses  and  they  will  be  denoted  by  er(i)  and  0(2).  The  lines 
which  are  perpendictilar  to  the  eigenvectors  of  this  matrix  aure  called  the  principal  lines. 

In  the  case  of  plane  stress,  the  third  principaJ  stress  0(3)  =  0  and  in  the  case  of 
plane  strain,  0(3)  =  i'(o(i)  +  0(2))«  Now  the  equivalent  stress  is  defined  as  follows: 

=  2  [(*^(1)  “  +  (^(3)  ~  ‘^(3))*  +  (^(3)  ”  <’■(1))^]  • 


5  Numerical  Results 


In  [28]  and  [30]  some  comparisons  were  made  between  MAM  (the  Method  of  AuxiliMy 
Mapping)  and  some  of  the  best  of  alternative  methods  such  as  Finite  Difference,  Finite 
Element,  and  singular  function  methods.  As  benchmarks,  elliptic  problems  having 
corner,  jump  boimdary  data,  or  interface  singularity  were  considered.  In  comparisons 
given  there,  it  was  shown  that  MAM  virtually  requires  no  extra  cost.  Since  CPU 
time  comparisons  for  elasticity  problems  between  MAM  and  alternative  approaches 
are  essentially  the  s£une  as  those  in  the  previous  papers,  we  only  compare  accuracy 
versus  DOF(Degree  of  Freedom)  between  MAM  and  the  conventional  approach  in  the 
framework  of  the  p-version  of  the  Finite  Element  Method. 

In  the  first  two  examples,  the  performance  of  the  Method  of  Auxiliary  Mapping  in 
the  framework  of  the  p-version  of  the  Finite  Element  Method  will  be  tested  using  the 
elasticity  problems  whose  true  solutions  are  known.  In  this  section,  all  computations 
are  the  plane  stress.  Recall  the  auxiliary  mapping  is  defined  by 

6)  =  (r^  cx}s{fi6),  sin(^0))  (42) 

The  number  /3  is  called  the  mapping  size  of  the  auxiliary  mapping.  Throughout  this 
section,  “  With  Map  ”  stands  for  the  results  obtsdned  by  applying  MAM  on  Mesh  I(the 
initial  coarse  mesh),  shown  in  Fig.  5.1.  “No  Map  ”  stands  for  the  results  obtained  by 
the  standard  Finite  Element  Method  on  Mesh  I  without  applying  MAM. 

Example  5.1.  Suppose  the  tractions  are  free  along  the  boundaries  a  =  ±ir  in  Fig. 
3.1.  Then  by  using  a  similar  argument  given  in  §3.1,  the  smallest  eigenvalue  is  A  =  0.5 
and  the  corresponding  stress  functions  are 

<7.  =  Ar(^-'){(2-Q(A-|-l))cos((A-l)0) 

-(A-l)cos((A-3)0)}  (43) 

(Ty  =  Art^-'>{(2-l-Q(A-l-l))cos((A-l)0) 

-f(A-l)cos((A-3)0)}  (44) 

Txy  =  A/^“^^{(A  —  l)sin((A  —  3)0) 

-|-g(A-bl)sin((A-l)0)},  (45) 

where  A  =  0.5  and  Q  =  1/3. 

Let  us  consider  the  equations  of  elasticity  on  a  domain  fli  =  {(x,y)  :  — 2  <  x  < 
2,  — 2  <  y  <  2}  shown  in  Fig.  5.1,  with  crack  along  the  negative  x-axis,  which  is 
isotropic  with  material  constants  E  =  1000(modulus  of  elasticity)  and  u  =  0.3  (Pois¬ 
son’s  ratio). 

The  traction  functions  given  by  (43)- (45)  are  imposed  along  the  entire  boundary 
of  ill,  including  both  sides  of  the  crack.  Furthermore,  the  following  constraints  are 
imposed:  the  displacement  vector  at  (0,0)  is  fixed  and  the  y-components  of  the  dis¬ 
placement  vector  at  (2, 0)  is  fixed. 
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On  the  mesh  shown  in  Fig.  5.1,  MAM  is  applied  with  mapping  size  =  4.  In 
order  to  show  the  effectiveness  of  MAM,  the  true  stress,  the  computed  stress  obtained 
by  using  MAM,  and  the  computed  stress  by  the  standard  Finite  Element  Method  are 
compared  one  another.  In  Fig.  5.2,  the  stress  Oy  on  [0,1]  x  [— 7r,7r]  are  compared.  In 
Fig.  5.3,  the  shear  stress  r^y  are  compared  on  the  same  subregion.  From  Fig.  5.2  and 
Fig.  5.3,  one  can  see  that  MAM  is  quite  effective  in  handling  crack  singularity.  In 
other  words,  one  can  not  see  any  difference  between  the  true  stress  and  the  computed 
stress  obtained  by  using  MAM.  However,  there  is  significant  difference  between  the 
true  stress  and  the  computed  stress  by  the  standard  FEM.  The  stress  tensor  along  the 
line  y  =  0.0129  are  given  in  Table  5.1. 

The  crack  singularity  discussed  in  Exaunple  5.1  is  not  too  strong.  Thus  it  is  possible 
to  obtain  a  practical  solution  by  sufficiently  refining  the  mesh  on  the  domain  Qi.  How¬ 
ever,  as  mentioned  in  §3.1,  there  are  some  elasticity  problems  containing  singularities 
which  are  too  strong  for  the  mesh  refinement  method  to  yield  any  practical  solutions. 

In  the  next  example,  it  will  be  shown  that  MAM  can  give  an  accurate  solution 
even  for  those  problems  to  which  the  mesh  refinement  method  alone  can  not  yield  any 
practical  solutions. 


(-2.3) 


mm 


Fig.  5.1.  The  scheme  of  the  domain  Ui  and  Mesh  I. 


Table  5.1.  The  Stress  ay,  r^y  along  the  line  y  =  0.0129:  True(True  stress),  With 
No  Map(computed  Stress  by  the  standard  FEM)  and  With  Map  (computed  Stress  by 
using  MAM).  The  .computed  stresses  are  obtained  by  using  basis  functions  of  order  8. 


<Tx 

X 

y 

True 

No  Map 

With  Map 

-0.40000 

0.0129 

0.05093822 

0.04823590 

0.05089369 

-0.34839 

0.0129 

0.06264135 

0.15461398 

0.06262561 

-0.29677 

0.0129 

0.07962216 

0.22911948 

0.07964541 

-0.24516 

0.0129 

0.10593000 

0.20079910 

0.10599321 

-0.19355 

0.0129 

0.15069789 

0.05516336 

0.15078407 

-0.14194 

0.0129 

0.23883382 

-0.05363119 

0.23890443 

-0.09032 

0.0129 

0.46349907 

0.37045664 

0.46355612 

-0.03871 

0.0129 

1.48632567 

2.41385925 

1.48644065 

-0.01290 

0.0129 

3.83449701 

4.62512498 

3.83435649 

0.01290 

0.0129 

4.42121220 

5.63512478 

4.42443179 

0.03871 

0.0129 

4.52331891 

4.53508576 

4.52399697 

0.09032 

0.0129 

3.25283283 

3.16439666 

3.25273651 

0.14194 

0.0129 

2.62989437 

2.47764699 

2.62995116 

0.19355 

0.0129 

2.26172620 

2.13848521 

2.26183925 

0.24516 

0.0129 

2.01336763 

1.95231425 

2.01344754 

0.0129 

1.57929020 

1.58108784 

1.57924127 

_ ^ _ 

Oxy 

True 

No  Map 

With  Map 

True 

No  Map 

With  Map 

0.00003312 

-0.02967877 

-0.00014424 

-0.00123211 

0.00510379 

-0.00123076 

0.00005368 

0.03406785 

-0.00003871 

-0.00173954 

-0.07977106 

-0.00176440 

0.00009401 

0.09624285 

0.00012811 

-0.00259535 

-0.13304302 

-0.00265020 

0.00018323 

0.08994510 

0.00034680 

-0.00417904 

-0.08602245 

-0.00426194 

0.00041800 

-0.03752552 

0.00063424 

-0.00752793 

0.08568925 

-0.00762896 

0.00123034 

-0.23020095 

0.00132026 

-0.01625616 

0.26892630 

-0.01635764 

0.00587308 

-0.15340532 

0.00583234 

-0.04945086 

0.07111119 

-0.04952390 

0.09964524 

1.07899882 

0.10064434 

-0.36327765 

-1.38006892 

-0.36333288 

1.83132606 

2.65316096 

1.83542437 

-2.41804125 

-3.04412791 

-2.41681990 

9.25729470 

6.04753022 

9.25774822 

1.00158548 

1.19976481 

1.00025011 

5.24987421 

4.78600541 

5.25097506 

0.69334021 

0.73212362 

0.69287164 

3.35173455 

3.25907239 

3.35178881 

0.22881300 

0.22695570 

0.22885594 

2.66240669 

2.52571398 

2.66266548 

0.11880174 

0.05174389 

0.11869870 

2.27678206 

2.17434706 

2.27711095 

0.07513994 

0.02115217 

0.07496224 

2.02172571 

1.97913221 

2.02198072 

0.05287338 

0.03560597 

0.05272972 

1.58175442 

1.57152579 

1.58166512 

0.02545255 

0.04008938 

0.02552800 
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Although  this  paper  is  concerned  with  the  equations  of  elasticity  on  polygonal 
domains,  for  brevity,  we  consider  an  elasticity  problem  on  a  sector  region  in  the  next 
example. 

Example  5.II.  Consider  the  equations  of  elasticity  (11)  and  (12)  in  a  wedge-shaped 
domain 

=  {(r,d) :  r  <  2, -a  <  ^  <  a},0  <  a  <  90° 

which  is  isotropic  with  material  constants  E  =  1000  and  u  =  0.3.  The  displacement 
functions  given  below  satisfy  the  equations  of  elasticity  in  the  domain  : 

Mr,d)  =  ^{-(A+l)/(0)}  (46) 

Mr, 6)  =  (47) 

where 

m 

A 

G 

The  corresponding  stress  tensor  are 

=  r("-'){(A-fl)-(A-|-l)2}sin(A-H)0, 

(Tg  =  r^^~^^(A  +  1)  sin(A  -|- 1)0, 

Trs  =  — A(A  +  1)  cos(A  -|- 1)0. 

It  is  worth  noting  that  ug  =  =  0  along  the  boundary  0  =  ±a.  Thus  this  problem 

could  be  solved  so  that  the  condition  =  0  would  be  imposed  by  proper  combination 
of  u*  and  Uy.  Nevertheless,  here  we  will  impose  the  displacements(nonhomogeneous 
essential  boimdary  conditions)  along  the  entire  boundary  of  the  domain.  Let  us  note 
that  the  displacements  are  strongly  singular  at  the  vertex.  In  this  example  we  show 
that  our  method  performs  very  well  even  for  this  case.  In  the  p-version  of  Finite 
Element  Method,  used  for  this  example,  the  essential  boundary  conditions  are  imposed 
by  the  L2  projection.  Hence  the  error  of  the  method  now  includes  also  the  error  of 
the  approximation  of  the  boundary  condition.  Therefore,  we  have  to  expect  that  the 
strain  energy  will  not  necessarily  be  monotonically  decreasing  with  p,  the  order  of  base 
functions,  as  it  would  be  if  the  nonhomogeneous  boundary  condition  would  be  satisfied 
exactly  (  See  also  remark  5.1).  Let  us  note  that  in  all  other  examples  mentioned  in  this 
paper  natural  botmdary  conditions  are  used;  hence  the  strain  energy  is  monotonically 
increasing  with  p. 


=  sin(A  -f- 1)0 


2(1  + 1/)' 
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Let  1^2“^  be  the  upper  half  of  as  shown  in  Fig  5.4.  Along  the  entire  boundary 
of  112**^  following  displacement  functions  in  the  x-y  coordinate  system  are  imposed: 

Ux  =  Ur  cos  6  —  U0  sin  9 
Up  =  Ur  sin  6  +  ue  cos  9 

which  correspond  to  the  displacement  functions  (46)  and  (47)  in  the  polar  coordinate 
system.  The  singularity  of  this  problem  can  be  as  strong  as  desired  by  taking  a  close 
to  90*  as  shown  in  Fig.  3.2  and  Table  5.2. 

Furthermore,  it  is  not  difficult  to  show  that 


=  AAr^^"^^sin(A  -  l)d, 
=  AAr^^“^^cos(A  —  1)0, 
=  AAr^^“')cos(A  —  1)0, 

OX 

dvx 
dy 


=  -AAr(^-*)sin(A-l)0, 


where  A  =  ~(A  +  1)/(2G).  Thus,  for  six  representative  wedge  angles  a  =  50®,  60®,  70*, 
80®,  85®,  89®,  the  true  strain  energy  (l/2)B({uej.},  {ue*}),  {ue*}  =  {u*,Uy}^,  on  the 
domain  are  as  follows: 


On 

On 

On 

On 

On 

On 


U{{u„))  =  0.4457011132335915D-02 
W({u„})  =  0.1531526418625024D-02 
W({ue*})  =  0.5573484518544861D-03 
W({u«:})  =  0.1707470731728283D-03 
U{{u„})  =  0.6899588499654107D-04 
W({uex})  =  0.1178312402203123D-04. 


Now  =  {(r,0)  :r<l,O<0<a}C  is  chosen  as  a  singular  region,  a 
neighborhood  of  the  singular  point  (0,0),  on  which  MAM  will  be  applied.  As  it  was 
mentioned  in  Remark  3.2, 

UgOip^  and  UyOifP  . 

are  imposed  along  the  boundaries  0*  =  0  and  0*  =  ^  of  (f2s“V  =  {(^*>^*)  •  ^ 

1,0*  <  a//9}.  Here  ^  is  an  optimal  mapping  size,  1/A  =  a/(90  —  a). 

Table  5.2.  The  eigenvalue  A  =  90®/a  —  1  for  six  representative  wedge  angles  a. 


a 

50® 

0 

0 

0 

0 

0 

0 

0 

00 

85® 

89® 

A 

0.50000 

0.286714 

0.125 

0.0588235 

0.01124 

Fig.  5.4.  The  scheme  of  Mesh  I  on  the  domain 


X-axis 


Table  5.3.  The  strain  energies  on  the  domain  for  a  =  50°,  60°,  70°,  80°,  85°,  89° 
when  MAM  is  applied. 


p-deg 

DOF 

50° 

60° 

o 

O 

1 

2 

0.438749274543 1457D-02 

0.1520859482894121D-02 

0.5567345729310525D-03 

2 

10 

0.4456887439527751  D-02 

0.1531572785147233D-02 

0.5573874898294548D-03 

3 

22 

0.445701 1922197191D  02 

0.1531527952481516D-02 

0.5573494680021682D-03 

4 

42 

0.445701 1143393066D-02 

0.15315264508271 15D-02 

0.5573484777837966D-03 

5 

70 

0.4457011132531040D-02 

0.1531526419381477D-02 

0.5573484525381 751 D-03 

6 

106 

0.445701 1132340058D-02 

0.1531526418643714D-02 

0.5573484518728477D-03 

150 

0.445701 1132336056D-02 

0.1531526418625500D-02 

0.5573484518549730D-03 

202 

0.445701 1132335964D-02 

0.1531526418625041D-02 

0.5573484518544919D-03 

p-deg 

DOF 

oo 

o 

o 

85° 

89° 

2 

0.1709561207603345D-03 

0.6906812121402299D-04 

0.11 7863953751 7020D-04 

10 

0.1707578614952583D-03 

0.6899855353039026D-04 

0.1178322837770610D-04 

22 

0.1707473741570929D-03 

0.6899596270247496D-04 

0.1178312715834583D-04 

42 

0.1707470816190254D-03 

0.6899588725323426D-04 

0.1178312411538886D-04 

5 

70 

0. 17074707341 16305D-03 

0.68995885061 98350D-04 

0.1178312402485034D-04 

6 

106 

0.1707470731796349D-03 

0.6899588499832774D-04 

0.1178312402218292D-04 

150 

0.1707470731730394D-03 

0.6899588499647640D-04 

0.1178312402211005D-04 

202 

0.1707470731728518D-03 

0.689958849964221  lD-04 

0.1178312402209949D-04 

Table  5.4.  The  strain  energy  on  the  domain  obtained  by  the  optimal  mesh 
refinement,  Mesh  V. 


p-deg 

DOF 

50° 

60° 

o 

O 

1 

10 

0.44056700205401 30D-02 

0.1566487115006998D-02 

0.6053103198646025D-03 

2 

42 

0.4457420260186373D-02 

0. 15344721 70765279D-02 

0.5647415490783892D-03 

3 

78 

0.4457064766350695D-02 

0.1531924078512398D-02 

0.5606895925910377D-03 

4 

138 

0.4457017168219616D-02 

0.1531597705602884D-02 

0.5597224836509342D-03 

5 

222 

0.4457011931428601D-02 

0.1531546547014034D-02 

0.5591 135803203382D-03 

6 

330 

0.445701124911 6450D-02 

0.1531535988017260D-02 

0.5587320835784172D-03 

.462 

0.445701 1150167903D-02 

0.1531531515488644D-02 

0.5582545129703971  D-03 

618 

0.4457011 13580271 lD-02 

0.15315294496331  lOD-02 

0.5579538617159367D-03 

p-deg 

DOF 

80° 

85° 

89° 

10 

0.2889017640162947D-)3 

0.3576618879140340D-03 

0.6182592096409667D-03 

42 

0.2534095506899352D-0;j 

0.3505035342285137D-03 

0.7129742782372550D-03 

78 

0.24731 1471 1669522D-03 

C.3607138018058794D-03 

0.7936052278053867D-03 

138 

0.2445714254542476D-03 

0.3701802501488918D-03 

0.8591464052345034D-03 

5 

222 

0.2337583505325639D-03 

0.337593145691 1586D-03 

0.7918101347523176D-03 

330 

0.2253173665961850D-03 

0.3105604155630747D-03 

0.7325557688296332D-03 

462 

0.2092777698103482D-03 

0.2447446052984088D-03 

0.5464419176391876D-03 

618 

0.1981310492965937D-03 

0.196880333421 1282D-03 

0.4054869521310821D-03 
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Throughout  this  paper,  to  measure  the  error  of  the  finite  element  solutions,  we  use 
the  following  definition: 


m\E.r  = 


W({«ex}) 


1/2 


(48) 


That  is,  it  is  the  square  root  of  the  difference  between  the  true  strain  energy  and  the 
computed  strain  energy  divided  by  the  true  strain  energy.  It  was  shown  in  ([36])  that 
||^||£,r  is  actually  the  relative  error  in  the  energy  norm,  provided  that  one  of 
the  following  cases  applies:  all  Boundary  conditions  are  either  homogeneous  Dirichlet 
or  arbitrary  Neumann  boundary  conditions;  some  Dirichlet  boundary  conditions  are 
nonhomogeneous,  but  all  other  boundary  conditions  are  either  homogeneous  Neumann 
or  homogeneous  Dirichlet  and  the  governing  equations  are  homogeneous.  Moreover 
all  examples  in  this  section  are  one  of  these  cases.  Hence  in  what  follows  we  will  call 
||^||js,r  the  relative  error  in  the  energy  norm. 

For  the  various  wedge  angles  a,  the  total  strain  energies  on  obtained  by  ap¬ 
plying  MAM  on  Mesh  I,  shown  in  Fig.  5.4,  are  given  in  Table  5.3.  The  relative  errors 
in  the  energy  norm(%)  are  given  in  Fig.  5.6.  By  comparing  with  the  true  solutions,  we 
can  conclude  that  MAM  is  able  to  yield  accurate  solutions  at  virtually  no  extra  cost, 
no  matter  how  strong  singularity  the  problem  contains. 

In  order  to  compare  the  results  obtained  by  MAM  with  those  obtained  by  the  mesh 
refinement  method,  Mesh  I  of  4  elements(see.  Fig.  5.4.)  are  refined  by  putting  circular 
layers  of  radii  <7,  <7^,  <r^,  <r^,  <r®,  <7®,  <t*,  centered  at  the  origin,  where  <7  =  0.15  (which 

is  known  to  be  an  optimal  geometric  ratio  for  a  geometric  mesh  refinenient  for  the  h-p 
version  of  FEM).  The  refined  mesh  obtained  by  putting  2,  3,  4,  ...,  9  layers  will  be 
denoted  by  Mesh  II,  Mesh  III,  Mesh  IV,...,  2md  Mesh  IX,  respectively.  These  meshes 
have  6,  8, 10,  12,...,,  20  elements  respectively.  For  example,  Mesh  V  is  as  shown  in  Fig. 
5.5. 

The  strain  energy,  obtained  by  applying  the  standard  FEM  with  the  refined  Mesh 
V  to  the  problems  on  the  domain  are  given  in  Table  5.4.  And  their  relative  errors 
in  the  energy  norm,  ||f  ||£7,r,  computed  by  using  the  true  strain  energy  are  shown  in 
Fig.  5.7. 

Remark  5.1.  The  energy  reported  in  the  tables  5.3,  5.4,  5.5  are  not  monotone 
with  p.  This  is  caused  by  the  fact  that  the  boundary  condition  is  imposed  only  ap¬ 
proximately.  For  higher  p  this  error  has  small  influence  and  the  energy  decreases  as  it 
would  occur  if  the  (essential)  boundary  condition  would  be  exactly  imposed. 

Table  5.4  and  Fig.  5.7  show  that  the  mesh  refinement  method  csin  handle  the  weaker 
singularities.  However  it  fails  to  give  any  practical  solution  of  elasticity  problems  with 
very  strong  singularities  such  as  the  cases  when  a  >  75°.  In  order  to  show  this  fact 
vividly,  the  standard  FEM  is  applied  to  the  case  when  a  =  89°  with  further  refined 


25 


meshes;  Mesh  I,  Mesh  III,  Mesh  V,  Mesh  VII,  and  Mesh  IX.  The  computed  strain 
energies  for  those  geometrically  refined  meshes  are  given  in  Table  5.4  and  Table  5.5 
and  their  relative  error  in  the  energy  norm  are  shown  in  Fig.  5.8.  In  this  case,  by  the 
massive  geometric  mesh  refinements,  finite  element  solutions  are  improved  a  little  bit. 
However,  the  accuracy  is  not  acceptable  at  all.  In  other  words,  the  h-p  version  of  the 
FEM  fails  for  this  problem. 

Remark  5.2.  (1)  Even  Though  ^  =  0/(90  —  o)  is  an  optimal  mapping  size 
for  the  problems  on  the  domain  fl2“\  MAM  with  other  choice  of  mapping  size  yields 
approximate  solutions  of  practical  accuracy.  Actually  MAM  with  mapping  size  (>  100) 
yields  the  relative  error  in  the  energy  norm  <  3%  when  the  order  of  basis  functions  is 
8  and  Mesh  I  is  used  for  the  problem  on  the  domain 

(2)  It  is  worth  noting  that  MAM  can  handle  the  elasticity  problems  even  when  the 
boundaries  of  the  neighborhoods  of  their  singularities  are  imposed  by  non-homogeneous 
essential  boundary  conditions. 


Fig.  5.6.  The  relative  error  in  the  energy  norm(%)  when  MAM  is  applied 
with  Mesh  I. 
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Table  5.5.  The  strain  energies  on  the  domain  '  obtained  without  applying 
MAM. 


p-deg 

DOF 

Mesh  I 

DOF 

Mesh  III 

2 

0.8348764761969085D-03 

6 

0.6720455342291412D-03 

10 

0.8363598689554478D-03 

26 

0.7753691417481764D-03 

22 

0.9184008854068470D-03 

50 

0.8631892906394158D-03 

42 

0.9931577497683493D-03 

90 

0.9345644101368223D-03 

5 

70 

0.9137775851488808D-03 

146 

0.8612360703348015D-03 

6 

106 

0.8430469500300520D-03 

218 

0.7967076857397165D-03 

150 

0.6223823176457848D-03 

406 

0.5940288770994173D-03 

202 

0.4550174531945867D-03 

510 

0.4405282879009032D-03 

p-deg 

DOF 

Mesh  VII 

DOF 

Mesh  IX 

14 

0.5688829276912590D-03 

16 

0.5235413474168798D-03 

58 

0.6556907603472391D-03 

72 

0.6030890129562007D-03 

106 

0.7297117726166816D-03 

132 

0.6710403386777638D-03 

186 

0.7898923442515647D-03 

232 

0.7262984181673099D-03 

5 

298 

0.7280587060296026D-03 

372 

0.66951 7682538331 7D-03 

6 

442 

0.6736470524085897D-03 

552 

0.6195529458233573D-03 

618 

0.5027442583422036D-03 

772 

0.4626180089851831D-03 

826 

0.3733095376025999D-03 

1032 

0.3437619833159763D-03 

NUMBER  OF  DEGREES  OF  FREEDOM 


Fig.  5.7.  The  relative  error  in  the  energy  norm(%)  on  the  domains  for 
a  =  50, 60, 70, 80, 85, 89  when  the  p-version  of  FEM  is  applied  with  Mesh  V. 
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10’ 


10* 


10* 


NUMBER  OF  DEGREES  OF  FREEDOM 

Fig.  5.8.  The  relative  error  in  the  energy  norm(%)  when  the  h-p  version  of  the 
Finite  Element  Method  is  applied  for  the  problem  on  the  wedge  domain  ^2“^  of  a  =  89®. 

In  the  first  two  examples,  the  smallest  eigenvalues  A  were  known.  Thus,  it  was 
possible  to  choose  an  optimal  mapping  size  =  1/A.  However,  in  engineering  practices, 
the  exact  eigenvalues,  which  represent  the  intensity  of  singiilarities,  are  not  known  in 
advance.  The  next  two  examples  demonstrate  that  MAM  succeeds  in  yielding  good 
approximate  solutions  even  for  these  cases. 

Example  5.III.  Let  us  consider  the  equations  of  elasticity  on  a  domain  H3  shown 
in  Fig.  5.9,  which  is  isotropic  with  materizd  const2mts;  E  =  1000  and  v  =  0.3.  Suppose 
the  boimdary  conditions  are  given  as  follows: 

(1)  u„  =  0,  «t  =  0  (fixed)  along  Fi  U  r2, 

(2)  Tn  =  10,  Tt  =  2  along  Fs, 

(3)  T„  =  0, 7(  =  0  (traction  free)  along  5n\(Fi  U  F2  U  F5), 

Then  it  follows  from  the  arguments  in  §3.1  that  this  problem  has  singularities  at  the 
crack  tip  Pi (0,0)  and  the  corner  P2(2,2).  Moreover,  suppose  the  body  force  is  zero, 
that  is,  {/}  =  {0,0}^,  then  min{Pe(Aj)}  are  approximately  0.3  and  0.7  at  Pi (0,0)  and 
P2(2,2)  respectively. 
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Fren 


Fig.  5.9.  The  domain  ils  and  the  line  segment  L  which  is  on  the  straight  line 
y  =  x(tan  ir/8). 

For  the  iinite  element  approach,  we  construct  two  meshes,  Mesh  I  and  Mesh  II 
on  n  as  shown  in  Fig.  5.10.  Mesh  II  is  obtained  from  Mesh  I  by  putting  layers  of 
radii  0.5,0.5<7,0.5<7^,0.5<t®  centered  at  Pi  and  layers  of  radii  0.5, 0.5(7,  0.5<t^  centered 
at  Pj,  where  <7  =  0.15.  Mesh  I  and  Mesh  II  have  22  and  48  elements  respectively. 
In  Example  5.III,  we  use  the  mapping  size  ~  %  and  =  2  on  the  singular  regions 
S\  =  {(ar,y)  :  ||(x,y)  -  (0,0)||  <  0.5}  and  S2  =  {(x,y)  :  ||(a:,y)  -  (2,2)||  <  0.5} 
respectively. 

In  Examples  5.III,  "  With  Map  ”  stands  for  the  results  obtained  by  applying  MAM 
on  Mesh  I.  **  No  Map  ”  and  ”  48EL  ”  stand  for  the  results  obtained  by  the  standard 
Finite  Element  Method  on  Mesh  I  and  Mesh  II  respectively  without  applying  MAM. 

The  total  strain  energy  obtained  by  the  following  three  ways  is  listed  in  Table  5.6: 
the  p-version  of  the  FEM  on  Mesh  I  by  applying  MAM(With  Map);  the  p-version  of 
the  FEM  on  Mesh  I  with  no  Map(No  Map);  the  p-version  of  the  FEM  on  Mesh  II  with 
no  Map(48EL). 
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Fig.  5.10.  (a)  Mesh  I  (22  elements);  (b)  Mesh  II  (48elements), 


Table  5.6.  Total  strain  energy  on  fta.  when  the  body  force  is  zero. 


p-deg 

With  Map 

48EL 

No  Map 

1 

1.702132029810480 

1.656404349217337 

1.231603737683330 

2 

2.061010428583605 

2.018717728104414 

1.605510783438384 

3 

2.093985634086386 

2.063320794501468 

1.736330296578600 

4 

2.110046396546211 

2.082748318204508 

1.815186503511321 

5 

2.113129774277219 

2.089615423359967 

1.864716288927604 

6 

2.113560741772014 

2.093441297022001 

1.899152183780876 

7 

2.113727076747088 

2.096139989799332 

1.925085281651720 

8 

9 

2.113785840066801 

2.113804271632455 

2.098163679729239 

1.945329260774251 

By  applying  the  extrapolation  approach  given  in  Chapter  4  of  ([36])  to  the  second 
column  of  Table  5.6  and  the  fourth  column,  DOF,  of  Table  5.7,  we  obtain  W(ttex)  = 
2.113815563245032,  the  computed  true  total  energy. 

Table  5.7  is  the  relative  error  in  strain  energy(%)  computed  by  applying  the  com¬ 
puted  true  energy  to  Table  5.6.  And  they  are  plotted  in  Fig.  5.11. 
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Table  5.7.  Relative  errors  in  the  energy  norm  (%)  and  Degree  of  Freedom  for  the 
cases  in  Table  5.6. 


p-deg 

With  Map 

DOF 

48EL 

DOF 

44.131 

64.603 

38 

46.517 

92 

15.805 

49.038 

120 

21.209 

280 

9.686 

42.259 

226 

15.456 

488 

4.223 

37.587 

376 

12.123 

792 

5 

1.801 

34.328 

570 

10.700 

1192 

6 

1.098 

31.867 

808 

9.818 

1688 

0.647 

29.880 

1090 

9.144 

2280 

0.375 

28.232 

1416 

8.605 

2968 

0.231 

1786 

In  this  example,  the  displacements  vary  within  very  small  ranges.  Hence  we  can  not 
see  a  clear  distinction  between  “With  Map”  and  “No  Map”  when  their  graphs  are  plot¬ 
ted.  However,  the  differences  are  clear  when  the  stresses  are  compared.  As  an  example, 
the  x-component  {U)  of  the  displacement  and  the  equivalent  stress  {a^)  at  the  points 
(r,  jr/8),  for  r  =  0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1,  and  -0.1,  -0.2,  -0.3,  -0.4,  -0.6,  -0.7 
along  the  line  y  —  xtan(7r/8)  are  listed  in  Table  5.8. 

The  equivalent  stresses  along  the  line  £(see.  Fig.  5.9),  given  in  Table  5.8,  for  three 
caiies  are  plotted  in  Fig.  5.12.  In  case  of  applying  MAM,  the  basis  functions  resemble 
the  true  solution  around  the  singularities,  hence  the  stresses  near  the  singularity  at 
/^(0, 0)  are  very  large.  Actually,  the  equivalent  stress  at  Fi(0, 0)  is  infinity.  We  plotted 
the  graph  of  the  equivalent  stresses  on  fl,  =  [-2,-2]  x  [—2, —0.1]  in  Fig.  5.13.  In 
Fig.  5.13,  the  x-grid  and  y-grid  sizes  are  0.2  and  0.1  respectively.  In  this  example,  the 
equivalent  stress  of  “With  Map”  is  much  bigger  than  “No  Map”  near  the  singularity 
at  Pi (0,0).  One  can  notice  this  fact  from  Fig.  5.13. 

The  equivalent  stresses  along  the  line  Z(see,  Fig.  5.9),  given  in  Table  5.7,  for  three 
cases  are  plotted  in  Fig.  5.12.  In  MAM,  the  basis  functions  resemble  the  true  solution 
around  the  singularities,  hence  the  stresses  near  the  singularity  at  Pi  (0,0)  are  very 
large.  Actually,  the  equivalent  stress  at  Pi  (0,0)  is  infinity.  We  plotted  the  graph  of  the 
equivalent  stresses  on  fl,  =  (—2,  —2]  x  [—2,  —0.1]  in  Fig.  5.13.  In  Fig.  5.13,  the  x-grid 
and  y-grid  sizes  are  0.2  and  0.1  respectively.  In  this  example,  the  equivalent  stress  of 
“With  Map”  is  much  bigger  than  “No  Map”  near  the  singularity  at  Pi  (0,0).  One  can 
see  this  fact  from  Fig.  5.13. 
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Fig.  5.11.  Relative  error  in  the  energy  norm(%)  for  the  domain  Q3  containing  two 
comer  singularities. 


Fig.  5.12.  The  equivalent  stress  along  the  line  L  shown  in  Fig.  5.9.  The  scale 
on  the  horizontal  axis  represents  the  radius  r  of  the  polar  coordinates  (r,7r/8)  of  the 
points  on  L.  For  these  computations,  basis  functions  of  order  8  are  used. 
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Fig.  5.13.  The  graph  of  equivalent  stress  on  fi.  =  [—2,2]  x  [—2, —0.1]  :  (a)  No 
Map:  (b)  With  Map.  The  x-grid  and  the  j/-grid  sizes  are  0.2  and  0.1  respectively.  The 
order  of  basis  functions  used  for  this  computation  is  8. 

Table  5.8.  Equivalent  stress  (<7*)  the  displacement  {U‘.  x-displacement )  along  the 
line  L  in  Fig.  5.9.  For  these  computations,  basis  functions  of  order  8  are  used. 


With  Map 


U 


6.836E-6  4.421 

2.516E-4  5.897 

5.075E-4  7.762 

7.695E-4  10.147 
1.289E-3  17.647 
1.535E-3  24.220 
0.2  1.754E-3  35.811 

0.1  1.908E-3  64.809 

-0.1  2.003E-2  44.053 
-0.2  2.253E-2  26.962 
2.430E-2  20.401 
-0.4  2.576E-2  16.809 
-0.6  2.831E-2  12.881 
2.949E-2  11.660 


No  Map 

U 

-2.153E-4 

3.945 

-1.243E-6 

5.234 

2.226E-4 

6.874 

4.501  E-4 

8.972 

8.953E-4 

15.504 

1.090E-3 

21.029 

1.250E-3 

31.129 

1.292E-3 

52.415 

1.547E-2 

36.843 

1.891E-2 

23.818 

2.072E-2 

20.083 

2.251E-2 

14.149 

2.524E-2 

12.135 

2.650E-2 

11.066 

48EL 


U 


-1.361E-5 

2.280E-4 

4.806E-4 

7.391E-4 

1.252E-3 

1.493E-3 

1.705E-3 

1.848E-3 

1.963E-2 

2.218E-2 

2.397E-2 

2.546E-2 

2.803E-2 

2.921E-2 


4.375 

5.833 

7.678 

10.038 

17.452 

23.936 

35.338 

63.761 

43.508 

26.716 

20.217 

16.687 

12.791 

11.636 
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So  far  we  have  considered  the  corner  singularities  on  isotropic  eiaistic  bodies.  How¬ 
ever,  the  interface  singularity  caused  by  an  abrupt  change  in  material  properties  in  zm 
elasticity  problem  has  a  similar  structure  to  that  of  the  corner  singularity.  However, 
the  interface  singularities  2u:e  usually  more  complex  and  stronger  than  the  corner  sin¬ 
gularities.  In  order  to  show  that  MAM  can  also  handle  this  type  of  singuljurity,  our 
next  example  concerns  an  elasticity  problem  with  two  interfaces. 


Example  5.IV.  Consider  the  equations  of  elasticity  on  the  domain  ^4  shown  in 
Fig.  5.14  that  is  composed  of  three  isotropic  materials.  That  is,  fl4  =  CI41  U  ft42  U 1143, 
and  the  material  constants  on  each  subdomain  are  given  in  the  following  table: 


On  041 

On  O42 

On  043 

E 

1000 

10 

1000 

1/ 

0.1 

0.001 

0.3 

We  also  assume  that  it  has  a  nonzero  body  force,  {/}  =  {10, 1000}^. 


Fix«d 


(2,2) 


Fr«e 


(2,0) 


Fig.  5.14.  The  domain  fl4  for  the  interface  problem  with  two  interfaces. 


This  interface  problem  has  singularities  at  Pi(0, 0),  ^2(2, 2),  and  P3(— 2, 2).  Thus  we 
choose  the  neighborhoods  of  the  singularities  as  follows:  Sj  =  {(a:,y) :  l|(x,j/)— P,||  < 
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0.5}, j  =  1,2,3.  A  mesh  on  U4  that  is  compatible  with  these  neighborhoods  of  the 
singularities  are  shown  in  Fig.  5.15. 

The  success  of  the  Method  of  Auxiliary  Mapping  depends  on  the  choice  of  the 
mapping  size  $  of  the  auxiliary  mapping.  In  our  method,  a  circular  sector,  S  = 
{(r,0) :  0  <  ^  <  <  To},  is  mapped  onto  S*  =  {(r*,5*)  :  0  <  <  9ol0,r*  <  r^^}, 

by  the  mapping  =  r^/^(cos0//9,sind/^).  Thus,  if  the  mapping  size,  is 

very  large  then  the  mapped  region  will  consist  of  very  narrow  circular  sector  elements. 
Nevertheless  the  convergence  theorem  given  in  the  previous  section  still  holds  since 
these  elements  satisfy  the  maximal  angle  condition  ([3])  which  allows  one  angle  to  be 
arbitrarily  small. 

In  this  example,  =  (A>^2>  A)  ”  means  the  results  obtained  by  applying  MAM 
vrith  the  auxiliary  mappings  of  size  for  the  singular  regions  81,82,83  respec¬ 

tively.  In  particular,  “/S  =  (1,1,1)”  stands  for  the  case  when  no  mapping  technique 
is  used.  The  total  strain  energy  obtained  by  the  various  choices  of  mapping  sizes  are 
listed  in  Table  5.9. 


Table  5.9.  Total  strain  energy  for  the  interface  problem. 


p-deg 

^  =  (8,10, 10) 

^  =  (5,10,10) 

DOF 

123032.2828452089 

123284.2678773375 

38 

150287.7980849437 

152583.4051500356 

116 

158804.9631362990 

159786.2535574225 

210 

161232.7426132401 

161314.8716586071 

344 

5 

161708.1813612756 

161697.6525341540 

518 

6 

161749.2509783793 

161745.8457528930 

732 

7 

161760.1348046330 

161758.1355019662 

986 

8 

161761.8532415435 

161760.3221783005 

1280 

9 

161762.9347649495 

161761.8908487467 

1614 

p-deg 

(9  =  (5,5,5) 

^=(1,1,1) 

DOF 

116027.58299108 

58106.6173116895 

38 

150578.57685495 

88079.9225407295 

116 

158821.08704690 

97828.5205550101 

210 

160167.37006738 

104140.5899691367 

344 

5 

160742.79149754 

108695.2364448717 

518 

6 

161051.55851408 

112211.9290644471 

732 

161239.18068546 

115041.7425766592 

986 

161361.96325978 

117386.4326246329 

1280 
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Fig.  5.15.  A  mesh  for  using  MAM  on  the  domain  CI4  with  two  interfaces. 


As  before,  extrapolating  the  second  column  of  Table  5.9  gives 

€iu„)  =  161765.7679239559 
as  the  computed  true  total  energy. 

Table  5.10  is  the  relative  error  in  the  energy  norm(%)  computed  by  using  this  true 
energy  and  Table  5.9.  And  these  relative  errors  in  the  energy  norm  versus  the  number 
of  degrees  of  freedom  are  plotted  in  Fig.  5.16. 

Table  5.10.  The  relative  error  in  the  energy  norm  (%)  and  Degree  of  Freedom  for 
the  cases  in  Table  5.8. 


p-deg 

=  (8,10, 10) 

^  =  (5,10,10) 

^  =  (5,5,5) 

DOF 

48.933 

48.773 

53.173 

80.050 

38 

26.637 

23.825 

26.300 

67.491 

116 

13.529 

11.062 

13.492 

62.869 

210 

5.740 

5.280 

9.940 

59.685 

344 

5 

1.887 

2.052 

7.952 

57.277 

518 

6 

1.010 

1.110 

6.645 

55.347 

732 

7 

0.590 

0.690 

5.705 

53.743 

986 

8 

0.492 

0.580 

4.996 

52.377 

1280 

9 

0.418 

0.490 

1614 
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Fig.  5.16.  Relative  error  in  the  energy  norm  (%)  for  the  interface  problem. 


Due  to  the  large  body  force,  we  can  see  a  big  distinction  between  ^  =  (8, 10, 10) 
and  =  (1, 1, 1)  even  in  the  displacements.  The  graphs  the  y-displacement  for  the  two 
cases,  =  (8, 10, 10)  and  =  (1, 1, 1),  are  plotted  in  Fig.  5.17.  Once  again,  the  stresses 
obtained  by  using  MAM  is  larger  than  the  stresses  obtuned  by  the  standard  FEM  at 
the  neighborhoods  of  singularities.  Numerical  experiments  show  that  the  singularities 
at  ^2(2, 2)  and  P3(— 2,2)  are  much  stronger  than  that  at  Pi(0,0). 

From  Table  5.9  and  Table  5.10,  we  can  conclude  that  the  mapping  size  ^  =  5  is 
not  big  enough  for  the  singularity  at  P2.  On  the  other  hand,  experiments  show  that 
the  mapping  sizes  =  10  and  ^  =  15  are  too  big  for  the  singularities  at  Pi  and  P2 
respectively. 


Remark  5.3.  In  all  cases  except  the  second  example  we  had  natural  boundary 
conditions  and  hence  the  energy  increases  with  p. 
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Fig.  5.17.  The  graph  of  the  y-displacement  over  fl  =  [—2,2]  x  [0,2]  :  (a)  No  Map; 
(b)  With  MAM.  The  i-grid  and  y-grid  sizes  are  0.2  and  0.1  respectively.  The  order  of 
basis  functions  used  for  this  computation  is  8. 


If  an  oversized  auxiliary  mapping  (i.e.  »  (min{i2e(A)})“^)  were  used  in  MAM, 

then  one  cannot  see  the  expected  improvement  until  the  polynomial  degree  is  of  high 
order.  It  W2is  shown  in  ([8])  that  if  larger  mapping  size  is  selected,  then  one  must  choose 
higher  degree  basis  polynomials  to  get  large  improvement  in  accuracy.  Hence,  from  a 
practical  point  of  view  an  optimal  choice  for  the  mapping  size  at  each  singular  region 
is  =  (min{i2e(Aj)})"'.  Thus,  in  order  to  obtain  optimal  results  from  MAM,  it  is 
desirable  to  know  the  eigenvalues  A  at  each  singularity.  Actually,  it  can  be  computed  by 
solving  trigonometric  equations  given  in  §2.1  for  the  corner  singularity  or  by  using  the 
computer  code  given  in  ([31])  for  the  interface  singularity.  Even  if  we  do  not  have  prior 
knowledge  of  the  eigenvalues,  MAM,  using  an  auxiliary  mapping  of  any  size  >  1, 
will  always  yield  a  large  improvement.  In  fact,  since  an  oversized  auxiliary  mapping 
yields  better  results  than  undersized  auxiliary  mapping  (unless  the  basis  functions  are 
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of  very  low  degree),  it  is  better  to  start  with  a  large  for  example  —  10  in  such  a 
case.  Another  possibility  is  to  use  different  strengths  of  mappings  and  select  the  one 
which  leads  to  the  largest  strain  energy. 


6  Concluding  Remarks 

MAM  can  efficiently  handle  the  plane  elasticity  problems  containing  such  singularities 
as  comer  and  interface  singularities.  No  matter  how  strong  singularities  the  problem 
contains,  MAM  yields  an  accurate  solution  at  vitually  no  extra  cost  if  the  str;> 'tures 
of  singularities  are  known.  Actually  MAM  can  handle  the  elasticity  problems  which 
even  can  not  be  solved  by  the  h-p  version  of  the  Finite  Element  Method.  In  applying 
MAM,  an  optimal  results  can  be  obtained  if  the  structure  of  the  singularities  are  known. 
However,  even  if  the  prior  knowledge  on  the  singularities  are  not  known,  MAM  can 
yield  very  reasonable  solutions  to  any  plane  elasticity  problems  with  singularity. 
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Tbc  Laboratoiy  for  Numerical  Analysis  is  an  integral  part  of  the  Institute  for  Physical 
Science  and  Tedmtdogy  of  the  University  of  Maryland,  under  the  general  administration  of  the 
Director,  Institute  for  Physical  Science  and  Technoi(^.  It  has  the  following  goals: 

To  conduct  research  in  the  mathematkal  theory  and  computational  implementation  of 
''  Clerical  analysis  and  related  topics,  with  emphasis  on  the  numerical  treatment  of 
ar  and  nonlinear  differential  equations  and  problems  in  linear  and  nonlinear  algebra. 

To  help  bridge  gaps  between  computational  directions  in  engineering,  physics,  etc.,  and 
those  in  the  mathematical  community. 

To  provide  a  limited  consulting  service  in  all  areas  of  numerical  mathematics  to  the 
University  as  a  whole,  and  also  to  government  agencies  and  industries  in  the  State  of 
Maryland  and  the  Washington  Metropolitan  area. 

To  assist  with  the  education  of  numerical  analysts,  especially  at  the  postdoctoral  level, 
in  conjunction  with  the  Interdisciplinary  Applied  Mathematics  Program  and  the 
programs  of  die  Mathematics  and  Computer  Sdence  Departments.  This  includes  active 
collaboration  with  government  agendes  such  as  the  National  Institute  of  Standards  and 
Tedmology. 

To  be  an  international  center  of  study  and  research  for  foreign  students  in  numerical 
mathematics  who  are  supported  Ity  foreign  governments  or  e»diange  agendes 
(Fulbright,  etc.). 

Further  information  may  be  obtained  from  Professor  L  BabuSka,  Qiairman,  Laboratory  for 
Numerical  Analysis,  Institute  for  Physical  Sdence  and  Technology,  University  of  Maryland,  CoUege 
Park,  Maryland  20742-2431. 


